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ABSTRACT 

This  report  investigates  the  feasibility  of  distance  determination  by  triangulation  for  artifi¬ 
cial  satellites,  rockets,  and  so  forth.  The  two-observer  case  is  solved  analytically  and  inves¬ 
tigated  in  detail  with  respect  to  the  propagation  of  errors.  For  objects  distant  compared  to 
the  observer-to-observer  distance  Ap,  the  variance  of  the  observcr-to-object  distance  R  is 
given  by  —  oR^  csc^  Ap.  a  is  the  standard  deviation  of  the  angular  measurements, 
and  0  is  the  angle  between  the  observer-to-observer  baseline  and  the  direction  to  the  object. 
The  next  topic  discussed  is  the  use  of  multiple,  simultaneous,  direction  determinations  for 
triangulation.  A  novel  method  is  proposed  to  deal  with  this  problem.  It  is  extended  to 
include  observations  with  different  measurement  precision,  and  generalized  to  take  into 
account  the  expanding  conical  nature  of  angular  errors.  Finally,  some  data  acquired  at  the 
Experimental  Test  System  of  the  GEODSS  network  and  the  Millstone  Hill  Radar  are 
analyzed  within  this  context.  It  is  clear  that  there  now  exists  a  powerful  observational  tool 
to  aid  in  initial  orbit  determination  utilizing  coordinated  angles-only  sensors. 
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DISTANCE  DETERMINATION  VIA  TRIANGULATION 


I.  INTRODUCTION 


The  problem  of  distance  determination  by  triangulation  arises  in  many  areas  —  surveying, 
passive  sonar  location  of  submarines,  parallax  determination  in  astronomy,  and  so  forth.  I  have 
examined  the  subject  anew  to  determine  its  utility  for  a  space-based  surveillance  system  which 
will  track  artificial  satellites,  rockets,  missiles,  and  so  on.  In  Section  II  of  this  report  the  ideal 
problem  is  solved,  t.ror  estimates  are  provided  for  the  results,  and  some  complications  of  the 
real  problem  introduced.  Section  III  discusses  a  novel  method  of  solving  the  problem  when  more 
than  two  observers  determine  lines  of  sight.  Although  designed  to  handle  multiple  observers,  a 
nice  feature  of  this  algorithm  is  that  it  reproduces  the  simple  solution  of  Section  II  when  there 
are  only  two  observers. 

These  two  early  treatments  assume  that  the  statistical  properties  of  the  measurement  errors 
of  the  different  observers  are  identical.  As  a  way  of  incorporating  varying  measurement  precision, 
a  pseudostatistical  version  of  the  problem  is  formulated  in  the  fourth  Section.  When  combined 
with  the  earlier  results  as  a  starting  point  and  an  iteration  procedure  designed  to  produce  rapid 
convergence,  this  final  technique  has  proved  to  be  robust,  easy  to  numerically  implement,  and  the 
best  of  several  methods  tried. 

The  last  Section  returns  to  the  two-observer  case  and  documents  the  analysis  of  data 
acquired  on  deep-space  artificial  satellites  at  the  Experimental  Test  System  of  the  GEODSS  net¬ 
work  and  the  Millstone  Hill  Radar.  It  is  clear  that  the  theoretical  expectations  have  been  veri¬ 
fied.  Hence,  a  powerful,  new  method  of  passive,  angles-only  localization  has  been  developed, 
tested,  and  polished  for  real  application. 


124826-N-01 


II.  THE  PURE  PROBLEM 


A.  FORMULATION 


An  object  in  space  has  a  geocentric  location  of  r  The  first  observer,  whose  geocentric  loca¬ 
tion  is  p|,  perceives  it  to  be  in  the  direction  of  the  unit  vector  Furthermore,  it  is  an  unknown 
distance  R|  from  this  observer.  It  follows  that  the  topocentric  location  of  the  object  is  R|  =  R|^| 
and  that 


Similarly, 


R|  =L-£i 


^2  -  ^2^2  -L~£2 


if  the  second  observer  measures  the  direction  to  the  object  at  the  same  instant  the  first  observer 
did*;  see  Figure  1. 


Figure  I.  Geometry  for  simultaneous  observation  of  an 
artificial  satellite  at  P  by  two  observers  at  O/  and  02- 
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*  Herein  lies  the  first  complication.  Because  light  travels  at  a  finite  speed  c,  even  if  both  observ¬ 
ers  are  carrying  synchronized  clocks  and  are  at  rest  relative  to  each  other,  and  if  they  both 
record  observations  at  times  t|  and  t2  and  furthermore  t]  =  t2.  they  still  have  not  made  simul¬ 
taneous  measurements.  Their  data  refer  to  locations  of  the  object  at  t]  -  R|/c  and  t2  R2  c, 
respectively.  Should  they  be  in  relative  motion,  the  whole  discussion  is  much  more  complicated. 


B.  SOLUTION 


With  reference  to  the  figure,  which  shows  the  two  observers  at  O]  and  O2  relative  to  the  geo¬ 
center  at  O  while  the  object  is  at  P,  let  <0|02P  =  02*  ‘^020|P  =  fli-  Then,  using  the  definition  of 

the  scalar  product, 

R2  •  (Pl  -P2>  -  ^2  1pi  -  P2I  COS  82 

Rl  ■  (P2  -  Pi)  =  Ri  IP2  -  Pll  cos  8, 

If  1  define £^2  the  unit  vector  from  the  first  observer  to  the  second,_pi2  "  (£1  “  P2)/IPl  ~  P2I’ 
then  1  can  rewrite  these  two  equations  as 

cos  8,  =  ^,2  •  «,! 

cos  82  =  p,2  -±2  (2) 

Since  1  knew  where  the  observers  were,  I  knew  what  directions  they  perceived  the  target  to  be  in, 
and  I  know  that  AO1PO2  ^i  plane  triangle  (at  the  instant  of  observation),  it  follows  that  both 
8|  and  82  can  be  uniquely  determined. 

I  can  now  determine  R|  and  R2.  From  the  law  of  sines  applied  to  AO1PO2  I  can  write 

R|  CSC  82  =  R2  CSC  8|  =  Ipi  -  P2I  CSC  [jt  -  (8]  +  82)] 

The  solution  for  R|  and  R2  is 

R|  =l£l  -£2lsin82  esc  (8,  +82) 

^2  =  !£i  -  P2ls‘n  ®i  CSC  (8,  +  82)  (3) 

Now  that  both  and  R2  are  known,  £may  be  computed  from  Equation  (1). 

C.  ANALYSIS  OF  VARIANCE 

Inevitably  there  are  measurement  errors  in  J,,  and  £2*  well  as  errors  in  £\  and  £2-  They 
can  be  propagated  into  values  for  the  variances  of  dy  and  82  and,  thence,  R|  and  R2  if  they  are 
small  and  random.  Because  of  the  I  — ’  2  symmetry  of  the  problem,  1  shall  only  discuss  the  var¬ 
iance  of  R|  =  var  (R|). 

There  are  three  groups  of  terms  in  the  expression  for  the  variance  of  By.  One  group  arises 
from  the  errors  in £y  and  £2  ‘heir  cross-correlations.  One  may  reasonably  expect  the  latter  to 
vanish  and  the  former  to  be  small,  equivalent  to  0701.  As  the  measurement  errors  associated  with 
^1  or_£2  will  be  at  least  an  order  of  magnitude  larger  than  this,  and  perhaps  a  thousand  times 
greater,  I  shall  ignore  all  errors  in  the  locations  of  the  observers. 

A  second  group  of  terms  in  the  expression  for  var  (By)  results  from  correlations  between  the 
observer’s  locations  and  the  measured  bearings  expressed  by  J,]  and  ^2-  While  systematic  errors  in 
the  observers’  locations  will  systematically  cause  and  ^2  incorrect,  there  is  no  reason  for 
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there  to  be  any  interplay  between  the  random  errors  in  the  observers’  locations  and  the  measured 
positions.  Hence,  all  these  terms  will  be  set  to  zero  also.  Of  the  36  elements  in  the  expression  for 
var  (0|),  only  three  are  left  (21  are  in  the  first  group  and  12  in  the  second). 

The  third  group  of  terms  has  its  origins  in  the  errors  of  measurement.  If,  to  be  concrete,  1 
write  as 

(cos  A I  cos  A  A 

cos  Ai  sin  Ai  j  ,  =  I  by  construction 

sin  A,  / 

where  A  is  a  latitudinal  angle  and  A  is  a  longitudinal  one,  then  the  reduced  expression  for 
var  (0j)  takes  the  form 


var  (0|) 


/  ad,  \2  /  36,  \3 


v.r(A,).2 


r)(^) 


cov  (A,.A|) 


The  covariance  of  Aj  and  A,  can  be  made  to  vanish  by  an  appropriate  observing  technique. 
Thus,  the  variance  of  6,  only  has  two  significant  terms. 

Utilizing  a  similar  chain  of  arguments,  one  deduces  that  the  dominant  terms  in  the  expres¬ 
sion  for  the  variance  of  R|  reduce  to 


var  (R,)  = 


/aRi\2  /3R,\2 


var  (62) 


since  cov  (0,,62)  =  0.  From  these  expressions  and  Equations  (2,3)  one  can  demonstrate  that 
var  (R|)  =  (£1  -  £2!^  ®2  csc^  (d|  +  62)  {cot^  (6,  +  ^2)  (^|) 

+  [cot  O2  -  cot  (6,  +  62)]^  var  (^2)}  (4) 

Let  us  try  and  clearly  see  what  this  means.  Interpret  A  and  A  as  the  usual  right  ascension 
and  declination  of  an  equatorial  coordinate  system.  Specialize  to  a  two-dimensional  situation  with 
Earth-based  observers  on  the  equator.  Finally,  let  p]  =  P2  "  P-  var  (A|)  =  var  (A2)  =  (for  sim¬ 
plicity),  and  set  AX  =  X2  -  >  0  e  -  the  longitude  difference  between  the  observers).  Now 

Equations  (2)  take  the  form 

2  cos  6,  -  [cos  (H|  -  AX)  -  cos  H|]  esc  (AX/2) 


2  cos  62  -  [cos  (H2  +  AX)  -  cos  H2]  esc  (AX/ 2) 


where  H  =  t  -  A  is  the  topocentric  hour  angle  (r  is  the  local  sidereal  time).  Furthermore. 


I  next  compute  the  variances  of  0j  and  ^2^ 

var  (0j)  =  (o/2)2  [sin  (Hj  -  AX)  -  sin  Hj]2  csc^  dj  csc^  (AX/2) 
var  (02)  =  (o/2)^  [sin  (H2  +  AX)  -  sin  H2]^  csc^  62  csc^  (AX/ 2) 

Finally, 

var  (Rj)  =  o2p2  sin^  62  csc^  (0,  +  62)  {csc^  0|  cot^  (0|  +  02) 

•  [sin  (H|  -  AX)  -  sin  Hi/p  +  csc^  $2  [cot  02  -  cot  (0)  +  02)]^ 

•  [sin  (H2  +  AX)  -  sin  H2f  } 

Even  this  expression  is  formidable.  As  a  last  simplification,  place  the  object  midway  between 
the  observers.  Then,  0|  =  02  =  0,  H]  =  -H2  =  H,  and  0  +  H  =  (rr  +  AX)/2.  Now, 

var  (01 )  =  var  (02)  = 

var  (R)  =  (o2p2/2)  sin^  0  sin^  (AX/2)  (csc^  0  +  sec^  0) 

This  is  still  not  transparent  but  becomes  more  so  as  the  object-to-observer  distance  greatly 
exceeds  the  observer-to-observer  distance.  For  as  r,  R  —  r  —  R  and  0  —  7r/2  from  below. 
Then,  R]  =  R2  =  R  =  2lpi  -  P2I  sec  0,  cos  0  =  (p/R)  sin  (AX/2).  Hence,  writing  0  as  nj2  -  60  for 
some  small  60, 

60  —  (p/  R)  sin  (AX/2)  ==  (p/r)  sin  (AX/2) 


var  (R)  - 


2p^  sin^  (AX/2)  2p2  sin^  (AX/2) 


As  a  numerical  example,  consider  two  GEODSS  sites  observing  a  near-stationary  artificial 
satellite.  With  AX  =  72°  (there  are  to  be  five  GEODSS  installations),  r  =  6.61p,  and  a  =  15",  the 
standard  deviation  of  the  topocentric  distance  would  be  17  km. 

D.  MORE  GENERAL  ANALYSIS  OF  VARIANCE 

A  generalization  of  Equation  (5)  to  arbitrary  geometries  would  be  useful.  1  provide  this  now 
only  for  the  case  when  the  object  is  far  enough  away  from  both  observers  that  the  sum  of  0|  and 
02  will  be  nearly  rr.  Indeed,  tt  -  (0]  +  02)  is  the  small  parameter  of  the  problem. 

Let^(a,6)  =  (cos  6  cos  a,  cos  6  sin  a,  sin  6)  be  the  unit  vector  of  direction  cosines  for  decli¬ 
nation  6  and  right  ascension  a.  As  before,  define  the  direction  from  observer  1  to  be  =^(a|,6i). 
The  pointing  vector  from  observer  2,  l2>  can  be  written  as 

I2  M 

for  some  vector  62,  of  small  norm.  Moreover,  because  22  is  a  unit  vector,  _2i  -62—0. 


From  AO1O2P  in  Figure  1,  the  relationships 

(£1  -£2)  'll 
|Pl  -  P2I 


cos  0 


cos  02  =  ■*■ 


(^1  -_P2)  '12 

1pi  -  P2I 


are  still  valid.  Setting  02  =  ”■  -  0]  -  60  these  expressions  take  the  form 


cos  0]  =  - 


(Pl  -£2)  'Al 
1pi  -  P2I 


with 


cos  02  =  -cos  (0|  +  60)  —  -cos  0|  +  60  sin  0| 


..  (P|-£2)ii 

60  =  — ; - ; CSC  01 


|P|  -£2! 

From  Equations  (3),  1  can  deduce  that  R2  is  given  by 
R2  =  [IPl  -_P2l  ®iP/[(^i  _P2)  iil 

and 

R|  =  R2  +  6R  ,  6R  =  -{pi  -_P2) 

As  the  variance  of  0|  is  simply  related  to  that  of  a|,6i  via 


var  (0,)  = 


301 
3a  1 


var  (0|)  + 


30, 


var  (6,) 


it  follows  that 

var  (0,)  =  [(p|2  •  3^,/3q|)2  var  (a,) 

+  (p|2  •  3SI|/36|)2  var  (6,)]  csc^  0, 


and  similarly  for  1  —  2  (where_p|2  =  (p,  -^2)/l£i  ~  £2!  still).  Examination  of  the  terms  in  the 
square  brackets  shows  that,  apart  from  the  var  (a,)  and  the  var  (6,)  terms,  the  remainder  is  some 
well-behaved  function  of  the  sines  and  cosines  of  a|,6]  and  the  direction  cosines  of^|2.  There¬ 
fore,  they  are  of  order  unity.  Assuming,  with  the  appropriate  metric  factor  subsumed,  that 
var  (a,)  ==  var  (6,)  =  Op  the  result 

var  (0|)  =  aj  csc^  0, 

follows.  Obviously,  this  is  just  as  true  for  the  variance  of  0^  to  terms  of  order  60. 
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The  last  formula  needed  is  Equation  (4).  Exploiting  the  fact  that  d|  +  62  =  ir  -  66,  and  that 
66  is  presumed  to  be  small,  the  leading  term  is  just 


i/i  o,R?csc2d, 

-K,  =  tv.-  (R,)]  -  I,, . 


(6) 


The  functional  dependence  of  0((|  on  R|,  |p|  -  pjl,  and  0|  are  the  key  elements  coupled  with  the 
fact  that  csc^  0|  ^  1.  The  trigonometric  terms  buried  in  Equation  (6)  cannot  play  a  major  role, 
nor  can  the  missing  \/Tin  the  expression  for  (2=1  +  1,  one  1  for  each  of  0|  and  02)- 


E.  COMPLICATIONS 


It  seems  unwise  to  assume  that  the  coordination  of  two  observers  will  be  so  successful  that 
they  will  be  able  to  make  truly  simultaneous  observations  (even  with  the  proviso  that  the  concept 
of  simultaneous  is  not  yet  precisely  defined).  A  much  more  realistic  scenario  is  the  following:  At 
some  time  t|  observer  1  commences  a  short  series  of  observations,  or  a  track,  on  the  object. 
These  independent  observations  are  separated  by  some  average  time  5t  and  are  Nj  in  number. 
The  last  one  is  performed  at  time  t2  =  t|  +  N|dt.  At  about  time  t|  (say  time  T|),  observer  2  com¬ 
mences  his  track  of  length  N26T,  terminating  at  time  T2  =  T]  +  N26T  which  is  near  time  t2.  We 
may  even  have  t2  <  T]  rather  than  t|  <  T]  <  T2  <  t2,  or  tj  <  T]  <  t2  <  T2,  and  so  on.  From 
these  two  tracks,  pseudo-observations  are  calculated  at  a  common  epoch  r.  This  computation 
involves  smoothing  the  two  tracks  in  some  statistical  fashion,  estimating  the  angular  velocity, 
the  angular  acceleration,  and  probably  higher  order  terms,  and  then  predicting  the  pseudo¬ 
observations  at  time  r.  This  will  be  a  substantially  better  procedure  for  short  observation  spans 
that  appreciably  overlap  than  for  short  observation  spans  that  do  not  overlap,  or  long  observa¬ 
tion  durations  in  general.  Therefore,  the  a  in  Equation  (6)  refers  not  to  the  sensor  pixel  size  but 
to  the  standard  deviation  of  the  pseudo-observation.  Very  roughly,  based  upon  the  simplest  least- 
squares  track  smoothing,  a  will  be  —  (pixel  sizc)l\/W if  the  observers’  locations  are  perfectly 
well-known,  and  so  on  (but  see  Subsection  V-A).  Otherwise,  one  should  include  the  hardware 
errors  and  discuss  an  "enhanced”  or  "effective”  pixel  size  which  is  larger  than  the  actual  pixel 
size. 


N  can  not  be  too  large,  for  then  N6t  will  be  too  long  and  the  track  smoothing  process  will 
be  limited  by  the  need  to  estimate  the  pseudoposition,  the  angular  velocity,  the  angular  accelera¬ 
tion,  and  so  forth.  As  the  VT* or  \/  3  or  even  \/T improvement  is  not  much  different  from  unity, 
there  may  not  be  much  gain  from  larger  values  of  N  at  a  fixed  dt. 

Should  the  observing  intervals  not  overlap,  or  not  even  be  nearly  centered  on  r,  then 
the  situation  rapidly  deteriorates.  The  sensitivity  of  the  variance  of  the  pseudopositions  to 
|t  -  (T|  +  T2)/2|  and  |r  -  (t|  +  t2)/2|  is  quadratic  at  first.  That  is  a  very  efficient  way  to  rapidly 
build  up  errors  that  will  swamp  (pixel  size)/\/^.  Thus,  while  my  “realistic”  scenario  relaxes  the 
perfect  synchronization  of  the  two  observers,  very  close  coordination  is  still  necessary.  I  shall 
continue  to  be  overly  optimistic  and  assume  this  to  be  the  case.  Thus,  a  —  (pixel  size)/ 2  (but  see 
Subsection  V-A). 
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Now  it  is  time  for  the  bad  news.  Look  at  Figure  2  which  shows  an  observer  making  two 
successive  measurements  of  a  moving  target.  Note  that,  because  the  observer  is  moving,  the  two 
directions  he  measures  are  in  different  coordinate  systems;  there  is  planetary  aberration  owing  to 
the  motion  of  the  observer.  Even  if  you  know  the  observer’s  motion  you  cannot  correct  for  this. 
You  must  also  know  the  object’s  distance.  The  majority  of  the  effect  is  wiped  out  (via  cancella¬ 
tion  of  opposing  terms)  if  the  epoch  of  smoothed  track  is  the  midpoint  of  the  observing  span. 
While  1  am  willing  to  be  optimistic  out  of  ignorance,  1  am  not  so  foolhardy  as  to  believe  that 
every  pair  of  tracks  will  have  the  same  midpoint.  Therefore,  in  addition  to  accidental  errors 
growing  as  |t  -  (t]  +  t2)/2p  and  |t  -  (T|  +  T2)/2|2,  there  will  be  uncorrected  systematic  errors 
growing  quadratically  with  these  time  differences  unless  a  very  complex  solid-geometry  problem 
can  be  solved.  So  far,  the  two-dimensional  version  of  it  has  resisted  all  attempts  to  solve  it. 

Is  this  a  fine  point  or  a  major  effect?  For  the  GEODSS  network  it  is  —  300"  (because  it  is 
composed  of  ground-based  telescopes  with  a  5t  of  2  min).  By  setting  r  equal  to  the  mean  epoch 
of  the  observations  all  but  2"  of  this  is  taken  care  of.  The  2"  residual  represents  the  unmodeled 
acceleration  of  the  telescope.  The  corresponding  amount  of  planetary  aberration  for  a  space- 
based  sensor  platform  could  easily  exceed  a  degree.  The  unmodeled  part  of  this  could  reach  Ol'5. 
As  long  as  instantaneous,  simultaneous  observations  represent  an  ideal,  this  source  of  error  is 
inherent  in  this  type  of  distance  determination. 
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Figure  2.  Illustration  of  planetary  aberration  owing 
to  the  motion  of  the  observer  (P/-"  P2)  and  the  object 
of  interest  (Q,  -  Q^J. 
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A.  FORMULATION 


S=(l/2)  1  li„-r„|2  (7) 

n,m=l 

n^m 

with  j;m  =  £fn  *  Rmim«  '*'***  vanish  for  error-free  correct  values  of  (Rm}-  In  practice, 

the  errors  in  the  observers'  locations  are  much  less  than  those  in  the  bearing  vectors  (see  above). 
Thus,  if  1  regard  S  as  a  “sum  of  squares,”  then  S  will  be  least  when  (R|,R2,...,Rn)  is  chosen 
subject  to  the  presence  of  errors  in 

B.  SOLUTION 

Since  S  is  a  quadratic  function  of  R,  if  I  form  the  “normal  equations”  by  V|^S  =_0,  then  1 
shall  have  a  system  of  N  linear  equations  to  solve.  In  particular,  this  system  is  just 

N  -  N  , 

(N  -  1)R„  •  S  R^  =in  •  (  -  N_p„  ,  n  =  1,2 . N  (8) 

m=l  V  m=l  / 

m^n 

While  the  above  is  cast  in  the  form  of  a  least-squares  discussion,  it  is  not.  There  is  no  sum 
of  squares  of  random  errors  being  minimized,  no  unbiased  estimator,  and  so  on.  Equation  (7) 
appeals  to  my  intuition,  leads  naturally  to  the  simple  Equation  (8),  and  will  be  shown  to  have 
some  geometrically  appealing  special  cases  (see  below).  Other  than  that,  its  justification  must  rest 
on  its  utility. 


III.  THE  MANY-OBSERVER  PROBLEM 


The  object  of  interest  is  at£.  There  are  n  =  I,2,...,N  observers  located  at  p|,P2.-”.Pn- 
Observer  N  measures  the  direction  to  the  object  and  expresses  the  results  as  a  unit  vector  of 
direction  cosines  The  “equations  of  condition”  arc 

_r=jOn  +  Rn|^  ,  n=l,2 . N 

for  the  N  unknown  distances  R|,R2,...,Ritj.  As  we  saw  earlier,  once  N  =  2  the  problem  is  soluble. 
Therefore,  for  N  >  3  the  problem  is  overdetermined  and  we  must  seek  some  statistical  fashion  to 
utilize  our  N  measures  of 

I  propose  to  do  this  by  exploiting  the  fact  that  the  object  can  only  occupy  one  place  at  one 
time.  In  particular,  the  quantity  S  defined  by 


C.  THE  TWO-OBSERVER  CASE 
If  N  =  2,  then  S  reduces  to 

S  =  RJ  +  r2  +  pj  +  p2  +  2R,i|  •  p,  +  2R2«:2  £2 
-  2(R|R2i|  •  ^2  *  R|ii|  ■  £2  *  R2i2  .£\  * £\  ’£2) 


From  dS/9R|  =  0  and  3S/9R2  =  0,  one  derives  (cos  6^^  'i.2) 

R]  -  R2  cos  0J2  =  -1pi  -_P2lP|2  ■  Ai 

-Rj  cos  fl|2  +  R2  =|Pi  -  P2IP12  ■ii2 

The  determinant  of  the  system  is  just  sin^  Hence,  there  cannot  be  a  solution  if  $12  =  0  or  tt 
(e.g.,  ^1  and  ^2  parallel  or  anti-parallel).  This  should  be  geometrically  obvious. 

The  solution  has  the  form 

Ri  =  -(ii  -I2  cos  012)  •^i2lPi  -JP2I  csc2  d,2 

R2  =  (^  -ii  cos  0,2)  -^pdfi  -  P2I  csc2  fl,2 

Note  that  J-i  -J,2  cos  0i2  is  perpendicular  to^,  ^2  -Ai  cos  0|2  is  orthogonal  toJ^j,  and  that  these 
are  merely  Equations  (3)  in  a  different  notation  [0|2  =  tt  -  0|  -  02  and  use  Equations  (2)].  The 
nonexistence  of  a  solution  in  the  intuitively  obvious  cases  and  the  reduction  to  the  deterministic 
solution  instill  confidence  in  the  meaningfulness  of  expression  (7). 

D.  THE  THREE-OBSERVER  CASE 

While  in  this  vein,  note  that  if  N  =  1,  then  Equation  (8)  just  reads  0  =  0.  More  interesting  is 
the  determinant  of  Equation  (8)  for  N  =  3.  With  cos  0„,„  =J,n  ■  it  may  be  written  as 

8  -  2  cos  0|2  cos  013  cos  023  -  2(cos2  fl|3  +  cos^  023  +  cos^  0i2) 

Once  again,  this  prohibits  a  solution  in  the  geometrically  clear  instances  —  two  0n,„  equal  to  tt, 
and  the  third  one  zero. 
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IV.  THE  FINAL  ALGORITHM 


A.  MOTIVATION 

Once  again  consider  an  observer  located  at  p  who  perceives  the  object  in  the  direction^. 

The  likelihood  function  for  the  object’s  location^is 

f(r)  =  C( wwO  exp  [-{r  -p)  ^Mfr  - jo)/  2]  (9) 

where  superscript  T  denotes  transpose  and  the  covariance  matrix  M  is  given  by 

M  =  wmm^  +  w'nnT 

C  is  a  normalization  constant,  while  w  and  w'  are  weights  that  provide  information  concerning 
the  precision  of  the  observations  along  the  directions  m  and  n.  These  two  unit  vectors  are 
orthogonal  toj..  If  the  measurement  errors  are  symmetric  in  the  plane  perpendicular  to^,  then 
M  has  the  form  of  a  projection  operator,  viz. 

M  =  w(mmT  +  nijT)  =  w(I  -g.g.T) 

The  likelihood  function  has  its  maximum  along  the  line  whose  tangent  isj^  and  is  constant 
on  all  parallel  lines.  Furthermore,  C  is  not  truly  a  normalization  constant,  for  Jf(r)d£  does  not 
exist  and  the  expected  value  of  is  not  defined.  If  !  adapt  this  model  to  the  two-observer  case 
then  these  difficulties  disappear.  Now,  f(r)  is  given  by 

f(jr)  =  Cf,(0  f2(r) 

where  each  likelihood  function  has  the  form  in  Equation  (9).  The  mean  value  for  x  is 
<X>  =  (M I  +  Mj)-'  (M  ,£,  +  M2P2) 

Its  covariance  matrix  is  just  (Mj  +  M2)'*.  This  matrix  exists  if  and  only  if  the  two  observers  plus 
object  are  not  collinear. 

For  N  >  2  observers,  the  process  is  simply  generalized.  Each  observer  has  his  own  location 
^n),  his  own  vector  of  direction  cosines  (£„),  covariance  matrix  (M^),  weights  (Wn,Wn),  and  so 
forth.  The  overall  likelihood  function  is  given  by  the  product  of  the  individual  likelihood  func¬ 
tions.  It  is  maximized  for 

j:=<r>=(2  Mn)  (  (10) 


whose  covariance  matrix  is 


Where  do  the  weights  w,w'  come  from?  Ordinarily,  weights  are  inversely  proportional  to  the 
squares  of  the  standard  deviations  of  the  measurement  errors.  In  this  problem,  that  would  not 
take  into  account  the  fact  that  a  fixed  angular  error  corresponds  to  different  linear  errors  at  dif¬ 
ferent  distances.  Hence  (see  Figure  3  for  motivation),  the  proposal  is  made  that  for  circularly 
symmetric  measurement  errors, 

w  =  w'=l/(Ra)2 

The  notation  is  as  above;  a  is  the  angular  measurement  error  of  an  observer,  and  R  is  the  topo- 
centric  distance  of  the  object  in  question.  Of  course,  the  R„  are  unknown  a  priori. 

BEARING  ERROR  FOLLOWS  N  (0.  TRUE 

=  sin  ’  (D/R)  DIRECTION 


OBSERVER  MEASURED 

DIRECTION 


Figure  3.  Motivation  for  w  =  i/(Ra)^.  If  the  bearing  errors  are  normally  distributed  with 
zero  mean  and  standard  deviation  a,  i.e..  as  N(0,a^),  then  the  displacement  errors  in  the 
tangential  direction  are  distributed  as  N[0.(Ra)^Jfor  small  a. 

B.  METHOD 

The  above  description  seems  a  reasonable  way  to  proceed  in  allowing  for  asymmetric  mea¬ 
surement  errors,  various  numbers  of  observers  with  different  measurement  precisions,  and  the 
incompleteness  of  the  angles-only  information.  However,  without  the  distance  information  the 
weights  can  not  be  computed.  Therefore,  I  recommend  the  following  procedure: 

(1)  Solve  for  r  by  using  the  method  of  Section  III.  Here,  r  is  understood  to  be 

N 

S  (R„i„  +  P„)/N 
n=l 

(2)  Use  the  set  of  topocentric  distances  and  the  known  measurement  errors  of  the 
observers  to  compute  the  set  of  weights  w„,w'„ ;  viz. 


W„=1/(R„o„)2  .  w'„=  l/(R„o'„)2 

If  the  measurement  errors  are  symmetric,  then  a„  =  a'„. 

(3)  Compute,  with  these  weights,  the  covariance  matrices  {M„}  and  then  the  solu¬ 
tion  for  X  from  Equation  (10). 

(4)  Calculate  new  values  of  the  topocentric  distances  from 

Rn=in(l'Pn) 

If  the  difference  between  these  topocentric  distances  and  the  old  set  is  suffi¬ 
ciently  small,  then  terminate.  Otherwise,  return  to  step  (2). 

This  method  is  straightforward,  numerically  simple,  robust,  has  an  obvious  convergence  cri¬ 
terion,  and  worked  better  than  several  other  algorithms  tested.  (The  details  of  the  numerical  test¬ 
ing  and  of  the  competing  algorithms  are  being  prepared  for  later  publication  elsewhere.) 
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V.  EXPERIMENTAL  TESTS 


The  Millstone  Hill  Radar  is  equipped  with  a  co-mounted  14-in.  telescope.  Having  a  wide 
fleld  of  view,  the  use  of  the  telescope  in  conjunction  with  the  acquisition  of  an  artificial  satellite 
speeds  the  radar  on  its  mission.  Unfortunately,  no  simple  means  exist  to  accurately  and  indepen¬ 
dently  calibrate  the  telescope.  Hence,  I  used  directional  data  acquired  by  the  radar,  rounded  to 
the  nearest  0.°01,  as  pseudo-optical  data.  This  precision  is  only  slightly  worse  than  that  obtain¬ 
able  at  the  ETS  in  its  Single  Star  Calibration  mode.  On  the  night  of  2  January  1984,  with  tele¬ 
phone  contact  between  the  ETS  and  the  Millstone  Hill  Radar,  a  short  series  of  tracks  were 
nearly  simultaneously  acquired  on  six  deep-space  satellites.  A  total  of  eight  sets  of  tracks  were 
obtained  with  four  to  nine  observations  per  track.  The  observing  frequency  was  once  every  30  to 
60  s.  These  data  were  smoothed  by  separate  least-squares  fits  to  second-order  polynomials  and  a 
pseudo-observation  generated  at  a  common  epoch.  These  pseudo-observations  were  and  £2  of 
Equations  (2).  The  common  epoch  was  the  mean  of  the  average  times  of  observation.  This 
can  be  shown  to  be  the  best  choice  when  the  two  sets  of  observations  are  equally  and  symmetri¬ 
cally  spaced  about  their  individual  midpoints.  Next,  knowing  the  locations  of  the  two  sensors, 
one  can  compute  P2,  P|2,  and  thence  d|  and  62-  Finally,  R|  and  R2  were  computed  from 
Equations  (3). 

Having  calculated  R  |  and  R2,  1  want  to  do  at  least  two  things:  assess  their  reliability  against 
the  “correct”  answer,  and  compare  the  magnitude  of  the  deviations  with  the  theoretical  analysis 
of  variance  presented  in  Section  Il-C.  What  is  the  “correct”  answer?  An  obvious  choice  for  R]  is 
the  distance  actually  measured  by  the  radar.  Such  a  comparison  is  shown  in  the  third  column  of 
Table  I  (the  units  are  kilometers).  Actually,  the  measured  distances  have  been  smoothed  with  a 
quadratic  polynomial  in  the  time  and  a  pseudodistance  generated  from  the  parameters  of  the  fit 
at  the  pseudo-observation  time.  Since  the  residuals  from  the  least-squares  fit  of  the  raw  distances 
are  generally  much  less  than  the  distance  mismatches  (see  the  last  column  of  Table  1),  this  artifice 
seems  harmless. 

A  second  choice  for  “correct”  might  be  the  distance  predicted  from  the  orbital  element  sets 
of  these  known  artificial  satellites  at  the  pseudo-observation  instant.  1  used  Space  Defense  Center 
element  sets  and  ETS  software  to  predict  Rj  values.  The  values  of  the  differences  between  these 
topocentric  distances  and  the  deduced  ones  are  in  the  fourth  column  of  Table  1.  There  is  a  very 
high  degree  of  correlation  between  the  values  in  the  third  and  fourth  columns  of  the  Table.  Since 
we  may  expect  that  the  Space  Defense  Center  gives  high  weight  to  Millstone  Hill  observations, 
and  that  the  ETS  software  is  correct,  the  correlation  between  these  values  is  not  a  surprising 
development.  Presumably,  had  1  utilized  Millstone  Hill  software  and  element  sets,  the  correspon¬ 
dence  would  have  been  even  tighter  (and  more  circular  in  its  logic). 

Can  one  predict  the  magnitude  of  these  discrepancies  from  the  analysis  of  variance  given  ear¬ 
lier?  In  general,  yes.  For  the  standard  deviation  of  an  individual  right  ascension  or  declination 
(actually,  horizon  system  coordinates  are  returned  by  the  radar)  one  can  not  use  the  canonically 
quoted  values  of  15"  for  the  ETS  and  18"  for  Millstone  Hill.  (However,  predictions  of  the  var¬ 
iance  of  Rj  utilizing  these  values  are  given  in  the  fifth  column  of  Table  1  for  comparison.  Note 


TABLE  I 


Values  of  and  Its  Errors* 


Satellite 

Track 


Ri 

Calculation 

RvMH 

Value 

RvEiSat 

Value 

15580.5 

71.1 

31.5 

25607.6 

-43.2 

15.2 

8817.8 

60.9 

19.5 

26591.1 

219.2 

188.8 

40761.4 

895.7 

901.9 

38090.9 

323.1 

310.0 

7659.5 

-21.9 

-71.8 

18364.6 

-659.7 

-577.7 

Standard  Deviations  of 


Theoretical 


Actual 


2243.5 

1515.4 
785.9 

3008.2 

3271.5 
10649.2 

1236.8 

1495.9 


All  quantities  are  in  kilometers. 


that  they  grossly  underestimate  the  actual  differences.)  A  better  value  for  the  standard  deviation 
of  an  altitude  a  or  azimuth  A  is  the  residual  from  the  least-squares  polynomial  smoothing.  Bet¬ 
ter,  but  not  best.  I  need  a  slight  statistical  diversion  to  further  elucidate  this  point.  However, 
before  delving  into  the  statistics  of  the  problem,  observe  that  the  sixth  column  of  Table  I  does 
contain  predicted  standard  deviations  utilizing  these  values.  Moreover,  the  magnitudes  are  gener¬ 
ally  in  consonance  with  the  third  or  fourth  columns  of  the  Table. 

A.  STATISTICAL  DIVERSION 

The  altitude  and  azimuth  (a,A)  generated  at  the  pseudo-observation  time  are  the  result  of  a 
statistical  adjustment  of  the  observed  altitudes  and  azimuths  {a„,A„}  and  their  times  of  observa¬ 
tions  {t„,t'„},  n  =  1,2,.. .,N.  This  explicitly  includes  the  computation  of  estimates  of  the  parameters 
Pt,P2,...,pvi  of  some  particular  model  fit  (a  polynomial  in  time  which  is  linear  in  the  p's).  A 
result  of  this  process  can  be  the  error  covariance  matrix  P  of  the  parameters. 

The  altitude  (or  azimuth)  at  the  pseudo-observation  time  can  be  written  as 

a  =  a(a|,a2,...,a^';  P|,p2 . Pm^  t|.t2.....t^;  fj.ts . t’^j 

The  predicted  Millstone  altitude  depends  not  only  on  the  observed  Millstone  variables,  the 


parameters  p  =  (P|,P2.— .Pm)»  observation  times  (t|,t2,...,tN),  but  also  on  the  times  of  obserxa- 

tion  of  the  same  satellite  at  the  ETS  (t'|,t'2 . t^)  for  the  pseudo-observation  instant  depends  on 

all  the  times.  Fortunately,  the  errors  associated  with  the  times  of  observation  are  so  small  1  can 
ignore  them  in  the  ensuing  discussion.  Thus,  all  times  are  considered  to  be  exact  and  the  explicit 
functional  dependence  of  a  or  A  on  them  will  not  be  exhibited. 

Now,  there  are  errors  of  observation  associated  with  the  ^a^,  An|.  Let  me  symbolize  this  by  a 
square  array  of  dimension  N  whose  main  diagonal  contains  the  variances  of  the  measurements. 
Call  it  V.  Therefore,  the  covariance  matrix  of  the  arguments  of  the  predicted  altitude  takes  the 
form 


The  covariance  matrix  of  a  and  A  has  the  form  of  a  product. 


wherein  (a|,a2 . a^;  Ai,A2,...,An)  is  a  vector  of  data  (in  practice  it  simplifies  matters  some¬ 

what  that  d  splits  into  two  independent  parts). 


If  1  write  out  the  above  matrix  products,  then  I  obtain 
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Examine,  for  example,  the  variance  of  the  azimuth  A.  It  takes  the  form 


2N 

=  s 

n=l 


(dAV  ,  /9A\t^/3A\ 
\9dn/  ‘^"^yap/  yap) 


The  first  term  represents  the  contribution  to  the  variance  of  A  from  the  random  errors  in  the 


measured  altitudes  (dj  =  a|,  d2  =  a2 . djM  =  aiyj)  and  azimuths  (d^+i  =  A|.  d]sj+2  =  A2 . d2\  - 


Aiyj).  The  curvilinear  nature  of  the  horizon  coordinate  system  implies  that  the  altitude  errors  have 
no  bearing  in  practice.  The  variance  of  the  n'^  element  of  d^  is  a^.  Also  note  that  each  element 
of  d  contributes  individually,  and  the  aggregrate  is  a  sum  of  squares. 

In  contrast,  look  at  the  second  term  in  a\.  It  involves  contributions  from  all  the  parameters 
of  the  fit.  Moreover,  the  derivatives  of  a  or  A  with  respect  to  p  may  themselves  be  functions  of  a 
or  A.  Thus,  the  “accidental”  errors  in  the  fitting  functions  create  systematic  errors  in  a  and  A. 
The  residuals  from  the  fit  contain  contributions  from  such  sources.  This  explains  why  the  values 


for  the  variance  of  R|  as  computed  using  these  values  are  so  much  larger  than  those  obtained 
using  the  “nominal”  position  uncertainties.  These  residuals  could  also  be  reduced  by  utilizing  a 
more  realistic  model  to  fit  the  data  to.  (The  answer  is  not  merely  a  higher-order  polynomial.  All 
the  computations  described  in  this  Section  were  also  performed  using  a  cubic  fit  —  with  marginal 
improvement.)  The  purpose  of  this  digression  is  to  explicitly  demonstrate  why  the  simple  “pixel 
size/x/N"  value  will  prove  to  be  a  significant  underestimate  of  the  true  angular  uncertainty. 
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